Theory of the collapsing axisymmetric cavity 
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We investigate the collapse of an axisymmetric cavity or bubble inside a fluid of small viscosity, like 
water. Any effects of the gas inside the cavity as well as of the fluid viscosity are neglected. Using a 
slender-body description, we show that the minimum radius of the cavity scales like ho oc t' a , where 
t' is the time from c ollapse. The exponent a very slowly approaches a universal value according 
to a = 1/2 + l/(4v/ — ln(t')). Thus, as observed in a number of recent experiments, the scaling 
can easily be interpreted as evidence of a single non-trivial scaling exponent. Our predictions are 
confirmed by numerical simulations. 

PACS numbers: Valid PACS appear here 



Over the last decade, there has been considerable 
progress in understanding the pinch-off of fluid drops, 
described by a set of universal scaling exponents, inde- 
pendent of the initial conditions 0, The driving is 
provided for by surface tension, the value of the expo- 
nents depend on the forces opposing it: inertia, viscosity, 
or combinations thereof. Bubble collapse appears to be 
a special case of an inviscid fluid drop breaking up inside 
another inviscid fluid, which is a well studied problem 
UQllijl: the minimum drop radius scales like ho oc t' 2 / 3 , 
where t' = to — t and to is the pinch-off time. Thus, huge 
excitement was caused by the results of recent experi- 
ments on the pinch-off of an air bubble Ll la l2l lill i or 
the collapse of a cavity ^lj m water, which resulted in 
a radically different picture, in agreement with two ear- 
lier studies 0, 0] . As demonstrated in detail in El ! 
the air-water system corresponds to an inner "fluid" of 
vanishing inertia, surrounded by an ideal fluid. 

Firstly, the scaling exponent a was found to be close to 
1/2, (typical values reported in the literature are 0.56 |9j 
and 0.57 [13), which means that breakup is much faster 
than in the fluid-fluid case, and surface tension must be- 
come irrelevant as a driving force. Secondly, the value 
of a appeared to depend subtly on the initial condition 
[Tl|. and was typically found to be larger than 1/2. This 
raised the possibility of an "anomalous" exponent, se- 
lected by a mechanism as yet unknown. To illustrate the 
qualitative appearance of the pinch-off of a bubble, in 
Fig. ^ we show a temporal sequence of profiles, using 
a full numerical simulation of the inviscid flow equations 
PJ. We confine ourselves to axisymmetric flow, which ex- 
perimentally is found to be preserved down to a scale of a 
micron [lfj, provided the experiment is aligned carefully 

a 

The only existing theoretical prediction 0, El El is 
based on treating the bubble as a (slightly perturbed) 
cylinder 0,0]. This leads to the exponent being 1/2 
with logarithmic corrections, a result which harks back to 
the 1940's llfj. Our numerics, to be reported below, are 
inconsistent with this result. Moreover, a cylinder is not 
a particularly good description of the actual profiles (cf. 
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FIG. 1: Numerical simulation of the time evolution of bubble 
pinch-off from initial conditions given by the shape with the 
largest waist. Pinch-off is initiated by surface tension, but the 
late stages are dominated by inertia, as observed experimen- 
tally [T3. 



Fig. 2J, as has been remarked before 0- In this Letter, 
we present a systematic expansion in the slenderness of 
the cavity, which is found to lead to a self-consistent de- 
scription of pinch-off. Our results are in excellent agree- 
ment with numerical simulations, and consistent with the 
experimentally observed exponents. 

Our approach is based on the standard description 
[171 Il8| of slender cavities, an assumption that is tested 
self-consistently by showing that the cavity's axial exten- 
sion is greater than its radius. The inviscid, irrotational, 
incompressible flow u = V(/> outside the cavity of length 
2L is written as 



(1) 



where C(£) is a line distribution of sources to be deter- 
mined. The length L will later drop out of the description 
of the pinch region, as indeed is not expected to be 
good near the ends of the bubble. For a slender geometry, 
d z <p <C d r <j), and the radial velocity, again using slender- 
ness, is easily evaluated to be d r <fi — —2C(z)/h(z). 

The equation of motion for the collapsing cavity of 
radius h{z, t) is dth w u r , and thus a{z, t) « — 4C, where 
a = h 2 and the dot denotes the time derivative. Finally, 
an equation of motion for C comes from the Bernoulli 
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equation, evaluated at the free surface |l£ 
arrive at 



We then 



fi(£,t)d£ 



a 
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4Ap/p, (2) 



where Ap = jk + const is the pressure difference across 
the cavity [2£|. In the two- fluid problem, the surface 
tension 7, multiplied by the mean curvature k ~ 
drives the problem. The capillary pressure will however 
turns out to be subdominant, so the last term in J5J) can 
effectively be neglected. Note that the resulting equation 
is invariant under a rescaling of both space and time, as 
both remaining terms are inertial (describing acceleration 
and convection of a fluid element). Thus dimensional 
arguments do not work, and a more detailed analysis is 
needed to fix the scaling exponent. Note that J5J does 
not conserve the volume of the cavity, whereas Fig. ^ 
assumes an incompressible gas inside the bubble. This 
however only affects the rounded ends of the bubble. 

Our aim is to explain the observed scaling behavior 
of the minimum cross section clq — a(0,i), as well as 
of the axial length scale A of the profile, which can be 
characterized by the inverse curvature A = (2ao/a' r ) 1 ^ 2 , 
where a ' = a"(0,t) and the prime denotes denotes a 
derivative with respect to z. Experiments as well as our 
own simulations show that ao as At l2a and A as Dt'° 
with (3 < a, thus the radius is small compared to the 
axial extend at the minimum. This means that a(0, t) 
can be neglected relative to £ 2 as A 2 in the denominator 
of the integral, except near the position £ = of the 
minimum. In other words, the integral is dominated by 
local contributions near the minimum. This will permit 
us to find equations of motion for the minimum in terms 
of local quantities alone. 

As shown later, a(£, t) goes to zero over the axial scale 
A. Thus the integral at z = can be approximated as 



/ [£ 2 + a r ^~a ln(2A 2 /a ). 



An arbitrary factor inside the logarithm depends on the 
exact shape of S(£, t); it can be determined empirically, 
but in fact becomes subdominant in the limit a ' — > 0. 
However, we now need another equation for the (time- 
dependent) width A to close the description. To that 
end we evaluate the second derivative of @ at z = 0. 
The contribution of the left hand side of (0) is 
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For a slender profile, a ' is subdominant, but the integral 
over the first term in angular brackets conspires to give 
zero in the limit ao — > 0, so the second term has to be 
considered as well, and a(£, t) has to be expanded beyond 
the constant term: a(£, t) — a + a '£ 2 /2. Thus using 
the same reasoning as before, and keeping in mind that 



0, we find for the second derivative of the integral 
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Equating this with the second derivative of the right 
hand side of J2J, (a 2 /(2a))", which is readily computed 
in terms of ao and A, yields the desired second equation. 
It is slightly more convenient to rewrite the results as 
equations for the time-dependent exponents 



2a = —d T ao/aa, 26 = — <9 T a '/a ', 



(3) 



where t = — hit' and = a — S. Note that is go- 
ing to be the "true" definition of the (time-dependent) 
exponents, which agrees with a local power-law fit. The 
result is 

(a r + a-2a 2 )ln(ri/a ') = -a 2 , (4) 
(S T + S - 26 2 ) ln(r 2 /a ') = 2a - 3a 2 - 2aS + 2a r ,(5) 

where the subscript denotes the r-derivative. 

The scaling factors ri,T2 have to be determined em- 
pirically, but only make a subdominant contribution as 
<Xq goes to zero. The time dependence of a ' is best found 
from integrating 



ln(a ') T = -2(5. 



(6) 



An analysis of Q-JBJ) shows that the approach to the 
singularity corresponds to an unstable fixed point as r — > 
00. As usual, this is the result of the freedom in the choice 
of singularity time to, see for example 0. The limit 
a = 1/2 thus has to be imposed onto the system in order 
to find the physically relevant solution. From the first 
equation, one finds that a approaches 1/2 from above, 
while the second equation says that 6 goes to zero, but 
remains positive. This guarantees the self-consistency of 
our approximation, although (3 approaches a in the limit. 
However, the approach of a and f3 toward their limiting 
values is exceedingly slow, as seen from the expansion 
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+ 0(r- 3 / 2 ), (7) 



where T is a constant which reflects the arbitrariness of 
the timescale in J2J) . Thus the value of T necessarily de- 
pends on initial conditions. However to leading order 
a approaches its limiting value in a universal fashion. 
Finally, for the self-consistency of our analysis we need 
that the dimensionless parameter a ' goes to zero toward 
pinch-off, as is indeed found from JBJ, owing to the slow- 
ness with which S converges toward zero. 

We now turn to a detailed comparison with full nu- 
merical simulations, not relying on any slenderness as- 
sumption, by focusing on the late stages of the pinch-off 



3 



event shown in Fig. To this end a suitably modified 
version of the boundary integral code developed to exam- 
ine inviscid droplet pinch-off was used, as originally 
reported in This involved two important modifica- 
tions: First, the boundary value operator (cf. Equation 
(If) in 0) has a zero eigenvalue in the case of the ab- 
sence of an inner fluid, corresponding to a change in the 
bubble volume. This singularity is analytically removed 
before the boundary integral operator is inverted, fixing 
the bubble volume. Second, due to the rapidity of bub- 
ble pinch-off, the adaptive time-stepping used for droplet 
pinch-off in |5j was replaced by a time-step halving pro- 
cedure with error estimation. 

A comparison of the numerical simulations with Q 
is given in Fig. Using equation ©, the value of a 
from the numerical simulations can be calculated as a — 
t'dfho/ho, and the pinch-off time to is estimated from the 
numerical data. The solid curve in Fig. 0is the data from 
the numerical simulation, the dashed curve is the leading 
order prediction given by equation Q with r = 0, and 
the dotted curve includes the adjustable constant with 

r = o.i. 

Data from the numerical simulations can be divided 
into three regimes. From approximately 10 -12 < t' < 
10~ 4 the bubble is considered to be in the asymptotic 
regime, and it is seen that there is very good agreement 
between the numerical data and the asymptotic theory: 
the leading order theory with T = accurately predicts 
the extremely slow decrease in the numerically deter- 
mined value of a, and the second order correction with 
r = 0.1 improves the agreement between the asymptotic 
theory and the numerical data. Equally good agreement 
was found for numerical runs using other initial condi- 
tions, provided that T was adjusted, as it is expected to 
depend on initial conditions, as also observed in experi- 
ment Time t 1 > 10~ 4 corresponds to a transitional 
regime where the bubble adjusts from an initial state 
where surface tension is required to initiate pinch-off, to 
an asymptotic state where surface tension is irrelevant. 
Time t' ~ 1CP 12 represents the threshold of the numerical 
simulations: extremely large interfacial velocities acting 
over ever- decreasing lengthscales, ultimately puts a limit 
on the validity of the numerical simulations. 

Gordillo et al. E3 have previously predicted that 
the minimum bubble radius h Q should scale with t' ac- 
cording to t 1 oc h 2 , y/— ln/iQ, using a method that in many 
respects is similar to ours |15J. However, the crucial dif- 
ference is that they do not treat the axial length scale A 
as a dynamical variable as we do, but effectively identify 
A with some outer length scale. Indeed, if one replaces 
a' ' by ao in (@J, one recovers the scaling result of [l5j . 
The conceptual difference between the two approaches 
is illustrated further by Fig. which shows the central 
peak of d from the full numerical simulation. The value 
of d rapidly drops to zero, effectively providing the cutoff 
of the integral J2J) at an axial length A, which is shrink- 
ing like t' 13 . So far, we have not been able to identify the 
logarithmic corrections of [3 in our full numerical sim- 
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FIG. 2: A comparison of the exponent a between full nu- 
merical simulations of bubble pinch-off (solid line) and the 
leading order asymptotic theory with T = (dashed line) and 
the second order correction with F — 0.1 (dotted line). 
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FIG. 3: A normalized graph of a = d 2 h 2 (z , t) / dt 2 as given by 
the full numerical simulations, for two different initial condi- 
tions, and at t' = 3.8 x 10" 10 (black line) and t' = 2.1 x 10~ 10 
(green line). 

ulations, since computing the axial scale is much more 
demanding than computing ho- 

In Fig. 0] we plotted the numerically computed min- 
imum radius ho, divided by the universal part of the 
present theory (full line), and that of 15] (dashed line). 
If normalized by an appropriate constant, the result 
should be unity. Namely, Q with T = is equiva- 
lent to h 0}Prei i cx t a / 2 \/ e~^~ lnt ' , while the theory in 
7] amounts to h 0iPred cx t a / 2 /{- In ft, 2 ,) 1 / 4 . While the 
present theory agrees extremely well with numerics with- 
out the use of any adjustable constant, the theory in [l5| 
varies by approximately ±50% over the range of t' plot- 
ted. 

In our earlier numerical simulations |l4j . as well as 
in most experimental papers 0, 0, ^| , the data for the 
minimum radius was represented by adjusting a single 
exponent a. Although Fig. [3 clearly shows that the ex- 
ponent is slowly varying, this subtle feature is difficult 
to detect in a more conventional plot like Fig. 0] To 
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FIG. 4: A normalized graph of ho/ho iPre d where feo.pred is 
predicted according to the theory presented by Gordillo et al. 
13 (dashed line), a least square approximation (dotted 
line), and the current asymptotic theory with a = 1/2 + 
(solid line). 



the data, as a result of the extremely slow variation of a. 
It also highlights the need for more sophisticated plots 
like Fig. [21 in the interpretation of future (experimental) 
data. 

To summarize, we have developed an asymptotic the- 
ory for the collapse of an axisymmetric cavity. A novel 
feature of this theory is a slow variation of the scaling ex- 
ponents, whose leading order contributions are universal. 
The slowness of the approach explains the experimental 
observation of apparently new scaling exponents, whose 
value may depend weakly on initial conditions. It re- 
mains to calculate the entire form of the central peak of 
a, which according to Fig. [31 is universal. This will de- 
termine the values of the constants Ti and T2- Other 
challenges are the inclusion of non-axisymmetry [t| and 
viscosity |10| into the theoretical description. 
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demonstrate this point, we have determined an effective 
exponent a — 0.559 from a least-square fit to the nu- 
merical data, a value which is close to those observed 
experimentally 0, ITo| . In essence, a can be viewed as 
the average over a values shown in Fig. O The resulting 
fit (dotted line) gives a surprisingly good description of 
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